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We compare difFerent quantum Master equations for the time evolution of the reduced density 
matrix. The widely applied secular approximation (rotating wave approximation) applied in com- 
bination with the Born-Markov approximation generates a Lindblad type master equation ensuring 
for completely positive and stable evolution and is typically well applicable for optical baths. For 
phonon baths however, the secular approximation is expected to be invalid. The usual Markovian 
master equation does not generally preserve positivity of the density matrix. As a solution we pro- 
pose a coarse-graining approach with a dynamically adapted coarse graining time scale. For some 
simple examples we demonstrate that this preserves the accuracy of the integro-differential Born 
equation. For large times we analytically show that the secular approximation master equation 
is recovered. The method can in principle be extended to systems with a dynamically changing 
system Hamiltonian, which is of special interest for adiabatic quantum computation. We give some 
numerical examples for the spin-boson model of cases where a spin system thermalizes rapidly, and 
other examples where thermalization is not reached. 



PACS numbers: 03.67.-a, 03.67.Lx 



I. INTRODUCTION 

The discovery that quantum computers would have 
much stronger capabihties for solving certain kinds of 
problems (such as number factoring (Ij or database search 
) than their classical counterparts has initiated a lot of 
research in quantum information theory 

Unfortunately, the fragile quantum coherence neces- 
sary for the superior performance of quantum computers 
is very sensitive to the inevitable interaction with the 
environment such that theoretical understanding of this 
process - called decoherence - is absolutely necessary [1] . 

For simple models (such as, e.g., a single spin or har- 
monic oscillator coupled to a thermalized bath of har- 
monic oscillators) and for sufficiently complex couplings 
to the reservoir after some time the system will equili- 
brate in a thermal state with the bath temperature 
This behavior would be consistent with our classical ex- 
pectations. 

A recent idea is to protect the quantum informa- 
tion by encoding the solution to a given problem in 
the (unknown) ground state of a problem Hamiltonian 
Hp . Given sufficient experimental control of the system 
Hamiltonian and a reservoir at sufficiently low temper- 
ature fcer <C ASniin (where AE^^i^ denotes the energy 
gap above the ground state energy of Hp), the ground 
state should be robust against decoherence in the sense 
that decoherence would always drive the system towards 
its ground state. In principle, one could then prepare the 
quantum system in any accessible state and wait suffi- 
ciently long until the equilibration has taken place. Un- 
fortunately, this mere-cooling approach is not expected 
to be very efficient, since the relaxation rate may be very 
small [H] or (in extreme cases) the system might get stuck 
in a local minimum [S]. 

A possible solution was proposed with the concept of 
adiabatic quantum computation 7]: The system is ini- 
tially subject to a simple Hamiltonian Hi and is prepared 



in its (easily accessible) ground state. Then, the Hamil- 
tonian is slowly deformed into the problem Hamiltonian 
Hp. The adiabatic theorem states that if this transfor- 
mation proceeds slowly enough, the quantum state will 
closely follow the instantaneous ground state Q. Fi- 
nally, for a nearly adiabatic evolution the system state 
approximates the system ground state to a high degree. 
Consequently, the maximum transformation rate (where 
the final excitations are acceptably small) corresponds to 
the computational complexity of the adiabatic algorithm. 
For closed systems, it is related to the spectral properties 
of the time-dependent system Hamiltonian M- For a 
reservoir at sufficiently low temperatures, this scheme is 
thought to be robust against decoherence Q and might 
even be aided by it fTlJ. 

Unfortunately, the standard framework of deriving 
master equations relies on some prerequisites that are 
not always fulfilled in realistic systems. For example, the 
Markovian approximation widely used is usually only for- 
mulated for time-independent system Hamiltonians. In 
addition, it sometimes leads to master equations that 
do not preserve positivity of the system density matrix 
[T^ [T3 |. Together with trace preservation, positivity 
grants stability of the density matrix eigenvalues and is 
thus necessary for probability interpretation (cf. [14]) 
of the density matrix. Consequently, observables ob- 
tained from non-positive density matrices may become 
unphysical 15]. This problem can be cured by the secular 
approximation. Combined with the Markovian approxi- 
mation, it leads to Lindblad-type [iB] master equations 
that generically preserve positivity of the density matrix. 
Unfortunately, the secular approximation is rather valid 
for quantum-optical but not for phonon baths Note 
that there exist non-Markovian master equations that are 
not of Lindblad type but nevertheless preserve positivity 
by construction. These models however are either phe- 
nomenologic [13, [3 in the sense that their parameters 
are not derived from a microscopic Hamiltonian or they 
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only grant positivity on a restricted set of initial states 
[19|. Especially in view of an experimental optimization 
of decoherence effects it is, however, necessary to relate 
the parameters in the master equation to those in the mi- 
croscopic Hamiltonian [20]. For example, for a realistic 
implementation of an adiabatic (or gate- model) quantum 
computer one would expect the qubits to be coupled to 
phonon degrees of freedom as well, such that a general 
treatment is advised. The present article shall present a 
further step in that direction. 

It has been noted [U, that coarse-graining may 

ensure for positive evolution of the reduced density ma- 
trix. However, the coarse-graining timescale so far had 
to be much larger than the inverse of the bath density of 
states cutoff. In this paper we argue that by adaptively 
changing the coarse-graining timescale one does not have 
to obey this constraint. Beyond this, we show that for 
infinitely large coarse-graining times we reproduce the 
widely used secular approximation. We show analyti- 
cally that for any fixed coarse-graining timescale t > 0, 
the resulting master equations are of Lindblad form, i.e. 
they preserve positivity and thereby also stability of the 
density matrix. 

This paper is organized as follows: In section [TT] we in- 
troduce our notation and in section |llT] we compare the 
standard procedure of deriving quantum master equa- 
tions from microscopic models with the proposed adap- 
tive coarse-graining scheme. We make our method ex- 
plicit by the example of the spin-boson model in section 
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II. GENERAL PREREQUISITES 



Denoting the time evolution operators of system and 
reservoir by U(t) and V(t), respectively, we can switch 
to the interaction picture 

Pit) = u\t)vHt)pit)v{t)uit), 

Hsjiit) = U\t)vHt)HsBV{t)U{t) 

A 

EE \Y,AMt)®B^it). (5) 

A 

We will denote all operators in the interaction picture 
by bold symbols throughout. In the interaction picture, 
the equation of motion for the density operator ([¥]) trans- 
forms into 

p(t)--z[i/sBW,p(i)] , (6) 

where one can exploit the smallncss of the coupling A to 
apply perturbation theory. 

Equation ^ can be formally integrated to yield 
t 

p{t) — p(0) — i J [H SB {t' ), p{t')] dt' and re-inserting this 



result in Eqn. ([6]) one obtains the following exact equa- 
tion 

p{t) = -z[/JsbW,p(0)] 
t 

- j [HMt),[HsB{t'),p{mdt' (7) 



for the density operator. 



We will consider Hamiltonians which can be divided 
into three parts 

H{t)^ Hs{t)+HsB + HB, (1) 

where Hs{t) describes the system part, Hb the part act- 
ing on the bath (with [Hs{t), Hb] — 0) and 

i?SB = A^A^®B^ (2) 

A 

denotes the interaction Hamiltonian with the small di- 
mensionless coupling parameter A <C 1 and system op- 
erators Aj, as well as bath operators Bj\ (differing cou- 
pling constants can be absorbed in the operator defini- 
tions). Hcrmiticity is only required for the complete sum 
(^^SB — -^sb)' splitting operators into hermitian 

and anti-hermitian parts one can always redefine them 
such that 

Aa = A\, Bj^ = b\, (3) 

which will be assumed further-on. 

The density matrix of the complete system is thought 
to evolve according to (fi = 1 throughout) the von- 
Neumann equation of motion 

p{t) ^-i[Hs (t) + HsB + Hb , p{t)] . (4) 



III. QUANTUM MASTER EQUATIONS 

We will first state the results of the Born-Markov ap- 
proximation without secular approximation (in subsec- 
tion |TTTX| and with the secular approximation (subsec- 
tion [1116]) in our notation. Afterwards, we will consider 
the coarse-graining approach in subsection IIII Cl 

With the usual assumptions (compare Appendix E]) 
involving initial factorization of the density matrix and 
neglecting any change in the reservoir part of the density 
matrix in ([7]) we obtain the Born equation 

PS = -*TrB{[-ffsB(t),Ps(0)p°]} 

t 

- I TTB{[HsBit),[HsB{t'),Psit')p°B]]}dt' 



+0{X'}, (8) 

where Tre {•} denotes the trace over the reservoir degrees 
of freedom. Evaluating the traces leads to the definition 
of the reservoir correlation functions 

C'ABiT) ^ TrB{e+'^--S^e-'^-^i?6PB} 

= C*BA{~r), (9) 
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and we obtain with {Bj\) = (which can always be 
achieved by a suitable transformation [2^ ) 



PS = / [[Asit')ps{t'),AAt)]CABit-t') 

+h.c.^ dt' + 0{X^} , (10) 



where h.c. denotes the hermitian conjugate. The integro- 
difFerential character of above equation complicates its 
solution, since analytical solutions are only possible in 
very simple cases [l^, HI] (see also subsection IIV Bp . and 
numerical solutions are hampered by the fact that the 
complete history of Ps{t') has to be stored in order to 
evolve ps{t)- 



A. Markovian approximation scheme 

In the usual Q Markovian approximation (see Ap- 
pendix |BJ one obtains for constant system Hamiltonians 
{Hs{t) — Hs) with the half-sided Fourier transforms 



(11) 



of the reservoir correlation functions ([9]) the time-local 
Born-Markov (BM) master equation (here given in the 
Schrodingcr picture) 



PS 



-i[Hs,psit)] 

+A2^^{r^g(^b-Sa) 



abed AB 

x{a\AB\b) {c\AA\dy 
+h.c.} , 



(|a) (6|)ps(i),(|c) (d\y 
(12) 



where Hg \a) — Ea\a) denote the orthonormal energy 
eigenbasis. By construction, Eqn. does preserve 

trace and hcrmiticity of pg. Note however, that positivity 
of its solution ps (t) is not generally preserved [13, [131 , see 
subsection IIV Dl for some counterexamples. 



B. The secular approximation 

In order to restore preservation of positivity in the BM 
approximation in general, it is necessary to perform the 
secular approximation (see Appendix [C|) . Typically, this 
approximation is known to be well-satisfied for quantum- 
optical systems, where it is also known as rotating wave 



approximation \2§\. In this case, one can combine 



+ 00 



IabH = TABi^) + nAi^)^ / C^B(r)e*""dr 



o-^s(w) = Tab((^) -Tba{^) 

CAB{T)sgn{T)e'^^dT 



(13) 



which yields the Born-Markov-Secular (BMS) approx- 
imation (in the Schrodinger picture) 

PS = -i[Hs,Ps{t)] 

^aab \a) {b\,ps{t) 

. ab 



+ EWrf[(|a> (6|)ps(<)(|c) {d\ 

abed 

-i|(|c) (d\)\\a) (&|),ps(i) 



'^ab.cd 



c AB 

x{c\AA\aY {c\AB\b) , 

AB 

x{a\AB\h) {c\AA\dY 



> jAB{Eb ~ Ea)5Ea-E^.Eb-Ea X 



(14) 



Eqn. (|14p has many favorable properties (compare e.g. 
chapter 3.3 in [3|): 

By construction, it preserves trace and hermiticity of 
the system density matrix ps. 

Since it is of Lindblad [i^| form (the matrix 7^b(w) is 
positive semidefinite), it preserves positivity of the den- 
sity matrix. 

For a thermalized reservoir characterized by the inverse 
temperature /3, the Fourier transforms of the bath cor- 
relation functions can be used to show that the system 
thermal equilibrium state with the same temperature 



Ps 



(15) 



is a stationary state. 

If the spectrum of the system Hamiltonian Hg is non- 
degenerate (implying that Se^^e^ — 5ab), the equations 
for the diagonal elements of ps in the eigenbasis of 
ffs completely decouple from the equations for the off- 
diagonals, and one obtains the same transition rates be- 
tween the populations as with Fermis Golden Rule. 



C. Coarse-graining approach 

Eqn. ^ is formally solved by 

p{t2) = 'W{t2,ti)p{ti)W\t2,ti) with the interaction 
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picture time evolution operator 



W{t2,h)^Texpl~i J HsB{t')dt'\ , 



(16) 



HsTi{ti)HsTi{t2)e{ti-t2) + HsTi{t2)HsTi{ti)Q{t2-ti), 

with Q(x) denoting the Hcavyside step function. Now, 
by expanding W{t + T,t) up to second order in A we 
obtain a second-order approximation to the fuU density 
matrix 



where the time-dependence of -H'sb(^) necessi- 
tates the time-ordering j28|] Ti?sB(^i)-f^SB(^2) = 



Pit + T) - Pit)~i 



t+T 



Hs-Biti)dh, p{t) 



t + T t+T 



J J {HsB{t2),Hsii{tl)]e{t2-h)dhdt2,p{t) 



t+T t+T 



Hss{tl)pit)Hs-B{t2) --{Hsii{t2)HsB{tl),pit)} 



dtidt2 + 0{X^} 



(17) 



Since such a truncated finite order approximation to the 
time evolution operator W(t + t) is still unitary, the 
above map (|17p preserves hermiticity, trace and positivity 
of the full density matrix p. An equivalent expression 
can be obtained from iterative solution of equation (O 
by keeping only terms up to 0{X^}. We will proceed 
with this second order approximation and derive from 
this fully unitary map a non-unitary map for the system 
part of the density matrix that preserves positivity. 

If we neglect the back-action of the system on the bath 
and assume factorization (Born-approximation) as de- 
scribed in Appendix [X] we can perform the partial trace 
over the reservoir degrees of freedom. By inserting the 
definition ^ in Eqn. (|17p and employing the definition 
of the reservoir correlation functions ([51 we obtain (again 
working in a frame where (i?^) — [22| ) 



Psit + T) 



Ps{t) 
A2 



t+T t+T 



O'Y. I I CAB{t2-h)sgn{t2 
X [A^it2)AB{h),Psit)]dtidt2 

t + T t + T 



h) X 



-ti X 



x^Ae{ti)ps{t)Aj^{t2) 
-l{Aj^{t2)As{ti),ps{t)} 



dtidt2 



Ps{t)+TL^{t)ps{t)+0{X^}., 



(18) 



which defines the action of the Liouvillian on the re- 
duced density matrix in the interaction picture. If 
one re-arranges the matrix elements of the density ma- 
trix as a A'^^-dimensional vector, the Liouvillian super- 
operator can be understood as an N'^ x N'^ (generally 



non-hermitian) matrix acting on ps- In the interaction 
picture the Liouvillian is small, i.e., since {Bj\} = we 
even have — 0{X^}. 

Unfortunately, the above map ps {t) Ps{t + t) has 
some shortcomings: It does not generally preserve posi- 
tivity of the reduced density matrix ps and in addition, 
if one applies above equation recursively n-times with 
small timesteps Ar such that nAr = r and compares 
the result with a single iteration of (HH]), the difference 
between the two solutions is much larger than OjA'^}, 
i.e., the solution depends on the choice of the stepsize. 

By defining time-averaging of an operator 0{t) over a 
time interval [i, i -)- r] as 



t + T 



\t,t+T] 



o 



we can write Eqn. p8|) as 



6{t')dt' , 



(19) 



// . \\ Ps{t + T) - psjt) 

{{PS))lt,t + T] = 

= L^{t)ps{t) + 0{\^}. (20) 

1. Explicit Liouvillian 

It is convenient to insert the even and odd Fourier 
transforms from ([13]) of the reservoir correlation func- 
tions 

-t-oo 

Cab(t) = ^ J lAB{uj)e-''^^duj, 

— OC 
+ 00 

C^6(T)sgn(T) = ^ J <7AB{uj)e-'^^duj, (21) 
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and to expand the system operators in the inter- 
action picture in the orthonormal energy eigenbasis 
Hs\a,) = Ea I a) of the system Hamiltonian (here assumed 
to be constant) 

ah 

Auih) = J2 \d) e^(^-^-)*i |c) {d\ . (22) 

cd 

I 



Then, we can use the relation 

t+T 

J e*"* dt' = re'"*e*"^/2sij^c [c^] ^^^h sine [x] = 
t 

together with Eqn. (PT|) and Eqn. in order to make 
the LiouviUian in Eqn. ((20)) exphcit. With denoting 
the energy differences of the system Hamiltonian by 
flab = Ea ~ E}y we arrive at 



Y,e^''-^'aab{T)\a) (6|,ps(i) 

, ab 

E e*^"=^"'^"'^*7c<i..a(r) [( \c) {d\)p^{t) ( l&> («l ) ^ - U ( \b) {a\)\ \c) {d\ ) , ps(t) 



abed 



47ri 



+00 



g.n„W2^^(c|^_^|„)*(c|^^|5) / ^_^g(,^)sinc[(u; + f7 



c AB 



{w + l^cfc): 



27r ' 



,i(0<:d-f26<,)T/2 



AB 



Y{b\AA\a)*{c\AB\d) / 7^e(^)sinc (w + f^ha) 



{w + fled); 



duj .(23) 



The coarse-grained derivative on the 1. h. s. of Eqn. 
(PO)) generates a time evolution that can be compared 
with its usual, first order differential equation counter- 
part (we denote the corresponding density operators by 
an overbar and the continuous index r referring to the 
LiouviUian chosen) 



dt 



(24) 



If we initialize this differential equation with a known 
density matrix pg and evaluate its formal solution at time 
t = T (i.e., after exactly the coarse-graining timescale), 
we find to first order in i"^ 



L^{t')dt' 



(25) 



This means that when considering the same initial con- 
dition, the difference between the two time evolutions 
resulting from Eqns. (|24p and ((20)) is given by 



Pi, (26) 



PlW-ps(r) = r^(r)+0{(L-)2}, 



7?(t) = 



T 

j L^{t')dt' -rL-'iO) 

.0 



i.e., essentially by the difference between the coarse- 
graining-averaged LiouviUian and its initial value. 

It follows from Eqn. ((^5)) that especially for short 
coarse graining times r ^ JaE[ — \ ('^ith Ai^max denot- 
ing the maximum energy difference of Hs) the difference 



will be negligible 77(t) « 0, such that the two solutions 
Ps{t) E^nd psiT) are equivalent. Since we have so far not 
made any assumption on separating timescales, this also 
implies that non-Markovian effects (that are expected for 
short times where the Markovian approximation does not 
hold) are within reach of the adaptive coarse-graining ap- 
proach, where the coarse-graining time is chosen to match 
the physical time. 

In the limit of very large coarse-graining times we also 
obtain lim ?7(t) « 0, since then the sine [•] functions in 

T — >00 

((^5)) begin to act like delta functions, see below. 

Also for intermediate coarse-graining times it is easy 
to see from the structure of the LiouviUian ((^5)) that the 
difference tj^t) is bounded throughout. 

From now on, we will omit the overbar and write Pg{t) 
for the solutions of Eqn. ((M)) . 



2. The infinite-T limit 

In the limit t — > oo one should note that for discrete 
a, b and under an integral over du with another integrable 
function one has in a distributive sense (see Appendix (F)) 



f{uj,a,b) — lim r sine {w + a) — 

r— »oc L 2. 

X 27r(5afc(5(w + a) , 



(27) 



where Sab is a Kronecker symbol and S{w + a) denotes 
the Dirac Delta distribution. Therefore, we obtain for 
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the Liouvillian matrix elements ([23|) in this hmit 
A2 



Icd.ba 



2i 



A^(5£;,-_E„,B,__E^ ^ \a)* X 



X (c| Ab \d)-iAB{Ea- Eb), 



(28) 



and we recover the secular approximation (I14p ! Natu- 
rally, this also implies that in the large-time limit, the 
solution Ps~*(i) captures all the favorable properties of 
Eqn. dUD. 



3. Positivity 

The most striking advantage of the coarse-graining 
procedure is that not only for very large, but for any 
fixed coarse- graining time r > 0, the resulting first-order 
differential equations are all of Lindblad form [l^ 
and thus intrinsically preserve the positivity of the den- 
sity matrix. This is easily seen by switching back to the 
Schrodinger picture, where the time-dependent phases in 
(l23t cancel 



^lcd.ba{T) 



abed 



(|c) {d\)pW){\b) {a\y 



l{{\b) H)^(|c) {d\),pUt) 



(29) 



which implies that the Schrodinger picture Liouvillian L"^ 
only depends on the coarse-graining timescale r. 

First of all, the first two commutator terms correspond 
to commutators with a hermitian operator, i.e., the sec- 
ond term accounts for the unitary action of decoherence 
(sometimes called Lamb-shift A]). Note that hermiticity 



of the corresponding effective Hamiltonian follows from 
the definition of the odd Fourier transform versus the 
half-sided Fourier transform (fT3|) . 

In order to show that p9|) is a Lindblad form, it re- 
mains to be shown that the matrix "fcd.baiT) is positive 
semidefinite. In order to see this we introduce the double 
indices i — (cd) and j — (6a) running from 1 to iV^ for an 
A^-dimensional system Hilbert space. Then, we can use 
the short-hand notation Wi — E^ — E^, Ag = (c| Ag 
Wj = Eh — Ea, A^ — (6 1 I a) and consider for arbi- 
trary complex-valued numbers xj in Eqn. ([23]) 



+ 00 



= ^ / ^z*y^{u;)-/ABiuj)zBi'^)duji30) 

where we have 

ZA = ^A^sinc[(w + Wj)T/2]e"^^^/2xj. (31) 

3 

Now, since the Fourier-transform matrix ^ab{'-^) is pos- 
itive semidefinite (compare also appendix [C| one has 
Zj({u!)jab{^)zb{^) > 0. The integral over a strictly pos- 
itive function can only yield positive results, and since 
also T > it follows that Jijir) — ^^cd.baij) is a positive 
semidefinite matrix. 

This also implies that solutions of the form (i) are 
always positive density matrices, since they correspond 
to an interpolation along Lindblad density matrices. 

D. Time-dependent Generalization 

Within the coarse-graining approach, it was in Eqn. 
([22|) where it was used for the first time that the system 
Hamiltonian Hs was time- independent and had a discrete 
spectrum. Here we will show that coarse-graining gener- 
ally leads to a time-inhomogeneous (i.e., one with time- 
dependent operators) Lindblad form of master equations 
and thus always preserves positivity of the density ma- 
trix. With introducing the notation 

l^(i,c^) = A^(i)e"^* (32) 

we can use equations (f2T|) and (fT9| in Eqn. (fT8|) to obtain 



((PS»[t,t+r] = 



47r« 



" 27r 



4-00 

I J E--m((^^h»|,,,.,((^-h»,,,/-'Ps(o 

-oo 

/ ^^lABi^^) 
AB 



Ab(^) X Ps{t)((AM 

t,t+T] \\ II t,t+T 



[t.t + T] 



(33) 
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With the replacements ((ps))[t j+t-j ^ p's"^ and 
Ps{t) ~> Ps(0 this becomes a time-inhomogeneous 
Lindblad master equation, since the positivity of the 
7^B(w)-matrix at every uj is guaranteed for reservoirs in 
thermal equilibrium. Intuitively, the time-dependence 
of the Lindblad operators does not destroy positivity, 
since at any fixed time t one may approximate the time- 
dependent operators by a constant-operator Lindblad 
form, see Appendix |E] for a more explicit discussion. 
Since the case of slowly varying system Hamiltonians is 
of special interest in the context of adiabatic quantum 
computation [^, we outline in Appendix |D] how one 
could in principle calculate the time-averaged operators 



t + T 



[t,t+r] 



t 

t+r 

- f U\t')AAUit')e+'^''dt' , 



(34) 



in the adiabatic limit. 



IV. SPIN-BOSON MODEL 



We will make our method explicit at the example of 
the spin-boson model in the following. We will also give 
some numerical solutions to the master equations used: 
The eigenvalues of the density matrices have been cal- 
culated with the LAPACK package The half-sided 
Fourier transforms ITTI) and odd Fourier transforms (fT5|) 
were calculated numerically from the full Fourier trans- 
form by consecutive application of backward and forward 
integral transforms. This was implemented by an inte- 
gration algorithm based on the discrete Fourier trans- 
form optimized for oscillating integrands (soj . For the 
integration of partial differential equations, a forth or- 
der Runge-Kutta method with an adaptive stepsize [s^l 
was used. Trace and hermiticity of the density matrix 
were always preserved within machine accuracy. Like- 
wise, whenever a Lindblad-type master equation was in- 
tegrated, non-negativity of the smallest eigenvalue of the 
density matrix was preserved within numerical accuracy 
defined by the accuracy of the Fourier transform. 



A. Microscopic Derivation 

We consider a system Hamiltonian with discrete energy 
eigenvalues that can be described by a quadratic form of 



Pauli matrices for a system of n spins 

n 



n n 



+EE E ^sv-f^ (35) 

i—l a.J3—x.y,z 

where ai denotes the Pauli matrices acting on spin 
i. In the worst case, this Hamiltonian is defined by 
l — 3/2n+9/2n^ real parameters. Note that when consid- 
ering explicit examples, we will not give the dimension 
of these parameters which implies that all times have 
their inverse dimensions. This system Hamiltonian is 
non-trivial in the sense that even in the case of time- 
independent parameters considered here, the time evolu- 
tion of the operators in the interaction picture cannot be 
solved analytically without using exponential resources. 
The system is coupled to a bosonic bath 



(36) 



with the usual bosonic commutation relations. The cou- 
pling between system and bath is realized by the quite 
general interaction Hamiltonian 



n 

HsB = E E 



Hk 



(37) 



i=l k 



with A <C 1 where the frequency-dependence is contained 
in the complex- valued coupling coefficients n^fc. 



1. Bath Correlation Functions 

In order to obtain a rather simple form of the master 
equation, we will make the assumption of a collective 
coupling, where the frequency dependence of the coupling 
strengths factorizes with the different spin positions and 
spin directions, i.e.. 



nik 



(38) 



for some function h^. This implies that the distance be- 
tween the spins is smaller than the correlation length of 
the reservoir oscillators. This approximation is not cru- 
cial for the further procedure but simplifies the resulting 
system of equations considerably, since the coupling to 
the reservoir can be described by just two effective spin 
operators. In this case, the interaction Hamiltonian (j37p 
can be written as 

i/sB = ASi^® [S + B^] +AS/®i [B-St] (39) 

with the composed operators 

n n 

Sfl ^5R(ni)-aj, Sj ^ (n,) ■ , 

i=l i=l 

B = J^hkbk, (40) 
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where we can identify from Eqn. ([2]) the hermitian oper- 
ators 



Bi ^ B + B\ B2^i(B- B^) 



(41) 



From the bath Hamihonian (|36p we obtain 
e'-f^BT^^g-j-ffBT ^ e"*"'=^6fe and the hermitian con- 
jugate, respectively. We will consider the limit of a large 
bath with a continuous spacing of oscillator frequencies. 
In this limit, the sums over k can be approximated by 
an integral 



J2\h,\'f{u;,) 



(42) 



and all other coefficients to vanish in Eqn. ([35)) as well 
as the simple coupling n = (0, 0, 1) in Eqn. ([ST]) - in this 
case also the collective coupling assumption (|38l) becomes 
exact. In this case the time evolution operator in the 
interaction picture (fT6)) can be calculated exactly 0, 
I33I I to all orders and one obtains that in the eigenbasis 
of the system Hamiltonian the diagonal elements of the 
density matrix remain unchanged and the off-diagonal 
elements simply decay as (cf. Eqn. (82) in [24[ in the 
limit of a continuous bath spectrum) 



T{t) = 



, ,sin2(wV2) , 
g{uj) \ coth 



duj , (47) 



where g{uj) is defined by the distribution of bath 
oscillators (spectral density) as well as the frequency- 
dependent coupling strengths hk- With the bosonic 
expectation value for a thermalized reservoir 

{Nk) = Tib {^I^fepJ^} = [eMM - 1]"' at tem- 

perature /? — (fceT)"^ this can be used to determine the 
bath correlation functions as 



i.e., in the pure dephasing limit one does not obtain ther- 
malization. Likewise, since g+^^st j^^-iHst = ^ it is ev- 
ident that BM and BMS approximations are equivalent 
for this case (compare also appendices IB] and IC)) . 



B. Non-Markovian Solutions 



Cii(r) 



Ci2(r) 



1 

2^ 
1 

2^ 



2^.9(1^1) 
|l_e-/3-'|' 



+ 00 



27r(-i)g(|w|)sgn(w)^_,^^^^ 



|l_e-/3^| 

C2i(t) = -Ci2{r), C22(r) = Cii(t), (43) 

where we can directly read off the Fourier transform ma- 
trix 



7(w) 



27rg(|u;|) / 1 sgn(tj) 
|1 - e-""^! V +isgn(tj) 1 



(44) 



Note that it is positive semidefinite (compare Appendix 
[C|) at every u ensuring positive evolution of the Lindblad 
form master Eqn. (|14p . 

In some examples given later-on, we will phenomeno- 
logically parameterize [sil, [s^l the density of states as 



(45) 



where the exponent s determines the behavior near small 
frequencies, Wph is some physical frequency of the bath 
and ujct is a cutoff-frequency necessary for normalization. 



2. Exact Solution for Pure Dephasing 

The limit of pure dephasing of a single qubit is defined 
by considering n = 1 with 



1 

7=2' ^ 



1 

'2 ' 



(46) 



For simple system and interaction Hamiltonians 



HsB = >^(y°- ® Bi 



where a € {x,y,z} and Bi = is a bath operator, 
the Non-Markovian Born Eqn. ^TU\i is analytically solv- 
able [H, [2^ [2^ in the special case of an exponentially 
decaying correlation function 

Cn(T) EE TrB {Bi(T)i3ipO } = ^e-l^l/^- . (49) 

By considering large Tb one can thus model reservoirs 
with a long-term memory, and the BM limit (where the 
BM-approximation becomes exact) is obtained by con- 
sidering lim Cii(t) = S{t). We obtain for the even and 

odd Fourier transforms in Eqn. (|13p 



(48) 



711 (w) 



1 + (^Tb)^ 



(Tll(w) 



1 + {LOTl,y 



(50) 



where it becomes visible (for later comparison with the 
spin-boson model) from Eqn. ((43|) that this case can in 
principle be reproduced by a bosonic reservoir in the large 
temperature limit with a Drude-like slowly decaying (but 
temperature-dependent) spectral coupling density 



(c.) 



2vr 1 + (cjTb)' 



(51) 



Note however that by assuming a sum of many exponen- 
tials and allowing for phases, the method can in principle 
be generalized also to the low-temperature limit [25l Is^ . 
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Inserting the operator definitions in the Born Eqn. ([10] 
we obtain 



and that the time evolution of the off-diagonal matrix 
element is governed by 



Vait) 



p-(t-t')/rb 

2^ 



dt' . 



PS = -[c7^ps{t)] + X^[rj,{t)-r^l{t),a^] , (52) 

where it becomes evident that by taking the time deriva- 
tive of the operator ry" one simply obtains a coupled set 
of first order differential equations 

PS = ^[a^ps(^)] + AM^a(^),^°] , 

Va = ^ [<J^fla{t)] - -fjait) + ^ K,ps(t)] ■ (53) 

for the operators f]a{t) = Va{t) — vHt) and ps{t)- Evi- 
dently, trace and hermiticity of ps and anti-hermiticity 
of fja are always preserved. Due to the initial condition 
Va{0) — 0, the trace of fja will always vanish. There- 
fore, it suffices to parameterize ps f^nd fj by just six real 
variables 



PS 



Poo Px + ipy 

Px - iPy 1 - Poo 

iVoo rjx + illy 
-Vx + iVy -T-Voo 



(54) 



In the limit Tb — > we simply obtain 



PS = 2 [a^ps(^)] + Wysith'^ - Psit)] (55) 

for the system density matrix. In the BM approximation 
(with finite Tb) we obtain along the lines of Appendix 
|B] for a = z (pure dephasing) rif^{t) = i(T^pg(i) and 
thereby 



d_ 



Pfit) \ 



-1/Tb 



-2A2 
-1/rb 



p'sHt) 



With the initial condition 77^^(0) = one finds 



pf{t) = p°i(0)e^*e-*/(2-^) coshfVl-SAVb-^ 

2Tb 



1 



, .- sinh -SAVb — 

VI - 8A^Tb V 2Tb 



t 



,(58) 



which reproduces the decay of the off-diagonal elements 
(the factor e** is a consequence of the Schrodinger pic- 
ture). In the high-temperature limit and for the corre- 
sponding density of states ([51]) leading to exponentially 
decaying correlation functions (^^ , the decay rate of the 
exact solution ([Tf|) becomes 



8A 



2 sm^{ujt/2) 



TT J Uj'^[l + (cJTb)' 




rduj 



4AVb 
— sinh 



cosh 
t 



(- 
I -J— I cosh , „ 

V2Tb/ V2Tb 



t 

2^ 



, (59) 



which reduces in the limit Tb — > to r(i) « 2A^i. Sim- 
ilarly, we find that in this limit, Eqn. ([58]) reduces to 
Ps^(^) = PgnO)e'*e~^'^ *. This can be understood as the 
limit Tb — > also corresponds to an infinitely fast relax- 
ation of the reservoir, where also the Born approximation 
becomes exact. Likewise, it is straightforward to see that 
\pfit)\ - e-r(*)pgi(0) = ©{A'*}. We will therefore not 
further discuss the pure dephasing case with exponen- 
tially decaying correlation functions and compare with 
the exact solution (|i7)) instead later-on. 



pgBM ^ ^ [a^^psit)] + X' [a^ps{t)a^ - Ps{t)] (56) 



2. Dissipative Coupling 



which coincides with (j55|) . i.e., the dependence 
on Tb vanishes. In contrast, for the more in- 
teresting dissipative case (a — x) we obtain 
V^'^it) = ^TT7^^"Ps(t) - ^T^^Vs(i), which yields 



Ps 



BM 



2 1+-r 



[[^^ps(^)],fT1 



A^ Tb 

'21+ri 



[[a\psit)],a^] 



(57) 



1. Pure Dephasing 

In the pure dephasing case one has a = z. Inserting 
the matrix elements in (|55|) one finds that pg° and are 
time- independent (as is also known from the full solution) 



Another important case is reproduced by choosing the 
dissipation coupling a = x in Eqn. (j52p . Inserting the 
ansatz ([5l)) into Eqn. ([55]) one finds two 3x3 systems 



/ Poo \ 








\Vyl 









2A2 




f poo 








Vx 





1 




{ Vy 




(60) 



which have an analytic solution that is too lengthy to be 
reproduced here. At first glance, these systems seem to 



be completely independent but note that the condition 
of initial validity of the density matrix relates their ini- 
tial conditions. The steady-state solution for the density 
matrix corresponds to the thermalized Gibbs state (flF 
for high temperatures (/3 ^ 0). 



C. Single-Qubit Coarse Graining 

1. Pure Dephasing 

In order to determine the master equations ((29|) for the 
pure dephasing limit discussed in subsection IIV A 2| we 
use 



fe 



\hkb. 



(61) 



to obtain that the Lamb shift Hamiltonian in Eqn. ((23|) 
is proportional to the identity matrix and thus has no 
effect. From the dissipative part however we obtain a 
non- vanishing contribution from Eqn. ()23p . such that 
the master Eqn. (|24p in the interaction picture is given 
by 

ps^ = f(r)[|0)(0|p-(<)|0)(0| + |l)(l|p-W|l)(l| 



|0) (0|p-(t)|l)(l|-|l) (1|p5W|0) (0| 



-f(T)p-(t) 

+ CXD 



f(r) 



27r 



7ii(a')sinc^ 



(62) 



With assuming exponentially decaying correlation func- 
tions (UH]) we can use Eqn. (|5D|) to recover the BM 
approximation (|56p . but above equation for pure de- 
phasing of course holds also for more general correla- 
tion functions. It is straightforward to show that the 
off-diagonal elements of the density matrix will decay as 
(0| plit) |1) = e-^(^)* (0| pliO) A closer inspection of 
the decay rate yields 



f(r)i = 8A2- 

T 



+ 00 



LU\ 

~2 |i_e-/5"| 



duj 



-t-CXD 



T 



-nr): 

T 



sin2(cjT/2) 



(63) 



i.e., we reproduce the result of [2J| that for pure dephas- 
ing of a single qubit, the solution p%^{t) yields the exact 
solution ps {t) , see also figure [TJ Naturally, this is also 
equivalent with the Born Equation up to 0{X^}. Note 
also that for large coarse graining times, figure [T] shows 
that the coarse-graining master equations approach the 
secular approximation master equation - which can also 
be seen as a numerical confirmation of Eqn. (I27p . 




10" 

time t [a.u." 



FIG. 1: [Color Online] Eigenvalue evolution of the sys- 
tem density matrix (only one eigenvalue shown since trace- 
conservation implies a symmetric evolution) initialized in the 
pure state (i| pg \j) = 1/2 for single-qubit dephasing. As pre- 
dicted by Eqn. (|63p . the solutions of the coarse-graining mas- 
ter equations (thin dashed lines) intersect the exact solution 
(solid line) at t = t (vertical dashed lines). The secular ap- 
proximation (thick dashed line) corresponds to r — > oo and 
does not correctly cover the short time behavior of the exact 
solution. Parameters were chosen as follows: (3 = 1, ujph = 1, 
ujct = 5, s = 1, X = 0.1. 



2. Dissipative Coupling 

With considering Ai = cr^, an exponentially decaying 
correlation function (|49|) . and Hs — ^[1 — a^] one ob- 
tains a more complicated master equation for p^{t) in 
the Schrodinger picture 



r' 



zaoo(r)lO) (0| - zan(T) |1) (1| , 



"701,01 (t) 



+7io,io(t) 



|0)(1|P^|1) (0|--{|l)(l|,p§} 
|1> (0|p^|0)(l|-i{|0) {0\,pl} 



+7oi,io(r) |0) (l|p^|0)(l| 
+710,01 (r)|l) (0|pg|l) (0| 



(64) 



From the non- vanishing matrix elements in Eqn. (|24p one 
can deduce that (unlike in the pure dephasing case) the 
Lamb-shift term will contribute, since (although diago- 
nal) it is not proportional to the identity matrix anymore. 
Likewise, we also observe here a decoupled evolution of 
diagonal 



Psoo(*) = 



7oi,oi 



1 



-(701,01+710,10)* 



701,01 + 710,10 
+e-('^oi.°i+^i«.i°)Vgoo(0) 



(65) 



and off-diagonal matrix elements 

Ps5i = - r(r)] pIq^ + 7oi,io(t) (Psoi)* , (66) 




time t [a.u.] time t [a.u.] 



FIG. 2: [Color Online] Evolution of the pn-matrix element 
of the density matrix for dissipative coupling. The solid line 
represents the solution of the Born Eqn. (|60p and the thin 
dashed lines correspond to solutions of (|64|) for different coarse 
graining times r. Up to second order in A, a perfect numerical 
agreement is found for the matrix elements of p'(f) (symbols) 
and the Born solution (solid line) . Parameters were chosen as 
follows: A = 0.1, Tb — 1. 

with r(r) = i (701,01 + 710,10) - i (o-ii - o-Qo) (we have 
omitted the r-dependence). Considering even and odd 
Fourier transforms of the form ([50|) corresponding to ex- 
ponentially decaying correlation functions, we can now 
compare the solution ((60)l of the Born equation with the 
solutions ([64]) of the coarse-graining approach. 

When one initializes the density matrix as 
Ps(0) = |0) (0| one finds that for infinite coarse- 
graining times the diagonal terms will equilibrate to 
value 1/2 (which corresponds to thermalization at 
infinite temperatures), see figure [D In this case, the 
off-diagonal terms will evidently vanish throughout. 
Note that for the small coupling chosen (A = 0.1), the 
solution of the Born equation ps(^) is approximated 
by the adaptive coarse-graining solutions pg(t) with 
extraordinary accuracy. 

In contrast, when initializing the density matrix as 
Ps(0) = i [|0) (0| + |0) (1| + |1) (0| + |1) (1|] one observes 
that the diagonal entries will remain unchanged and the 
off-diagonals will decay. The actual behavior of the de- 
cay of the off-diagonal elements is depicted in figure [S] 
Again we observe a strikingly good agreement between 
the Born solution psit) and the adaptive coarse-graining 
solutions Pg{t). 

It is also instructive to compare for the adaptive coarse- 
graining approach the limit of infinite coarse-graining 
times (BMS approximation) 

+ T^m{MPs\l) {0\ + \l) (OIpsIO) (1|-ps] 

(67) 



FIG. 3: [Color Online] Evolution of the |poi I matrix element 
of the density matrix for pure dissipation. The solid line rep- 
resents the solution of the Born Eqn. (|60p and the thin dashed 
lines correspond to solutions of (|64|l for different coarse grain- 
ing times T. Up to second order in A, a perfect agreement 
is found for the matrix elements of p^{t) (symbols) and the 
Born solution (solid line). The difference between Born so- 
lution and BM as well as BMS approximations is too small 
to be visible. Parameters were chosen as follows: A — 0.1, 
Tb = 1. 

with the BM approximation master equation (j57[) , where 
one can see that only the equations for the diagonals 
match. Likewise, in the limit Tb — > one obtains 

+AMIO) (l|pg|l) (0| + |1) (0|pg|0)(l|-pg] 
-t-A^e— sinc[T]|0)(l|p§|0)(l| 
-|-A2e+-sinc[r]|l)(0|p^|l)(0| . (68) 

Here, one finds that the BMS approximation (r — > 00) 
yields a different master equation than the BM limit 
(l55)l and also the BM approximation ((57| However, also 
within the coarse-graining approach the limit 
(with finite r) leads to the same steady state as the BMS 
approximation t 00 and the BM limit (|55p - only the 
relaxation rates may differ. Note however that the BM 
approximation (j57p may even lead to instable behavior 
of the off-diagonal matrix elements for large A (compare 
also section |IVD[) . whereas the coarse- graining approach 
always generates Lindblad forms. 

3. General Coupling 

For a more complex coupling of the spin to the reservoir 

Hs = i[l-a^] , ni = -i=(l + z,l + z,l + i),(69) 

and also more realistic spectral coupling densities of type 
we find thermalization as predicted by the BMS ap- 
proximation for long times. For short times however, the 
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time t [a.u.] 



A = 3?{r(+l) + r(-l)}. For physically motivated bath 
correlation functions Cn in Eqn. P5)) one obtains that 
(Ti < 0, such that the evolution in this subspace may not 
lead to unstable behavior - although positivity may be 
violated. 

In contrast, in the off-diagonal subspace one obtains 
the two eigenvalues A2/3 = -2}? A±^/ -\ + 2B\^ + A^X'' 
with A defined as above and B = 3{r(-l) - r(+l)}. 
Given i? > (which can be achieved with correlation 
functions of form ((431) '). oi^^ of these will pick a positive 
real part as soon as > (25)"^, which corresponds to 
an unstable evolution, see also figure El Numerically, we 
observe the same for n = 2 mutually uncoupled qubits 
(with similar system Hamiltonians and couplings). In 
this case, the expectation value of operators that are not 
diagonal in the system Hamiltonian will not yield any 



FIG. 4: [Color Online] Evolution of the eigenvalues of the 
density matrix for a single qubit initialized in the pure state 
(«|Psli> = 1/2. For long times the BMS solution (thick 
dashed line) and the adaptive coarse-graining solution (thick 
solid line) predict thermalization in the thermal equilibrium 
state (jlSp . whereas a different equilibrium state is assumed 
for short coarse-graining times (thin dashed lines). Parame- 
ters were chosen as follows: /3 = 1, cjph = 1, ciJct = 5, s = 1, 
A = 0.1. 

are chosen too small, a non-thermalized steady state is 
obtained. For large coarse-graining times, the thermal- 
ized steady state is reached, but then in the short time 
regime, large differences between fixed graining and the 
adaptive graining solution are found. Again, we numer- 
ically confirm Eqn. (|27p . since for large coarse graining 
times the secular approximation is reproduced. 



D. Staying Reasonable 

Frequently, the BM approximation is used al- 

though it does not generally guarantee for positive evo- 
lution. This is certainly tolerable if the solution obtained 
approximates the exact solution well (and thus violates 
positivity only slightly). In addition, if the interaction 
Hamiltonian used in such models implicitly corresponds 
to a secular approximation [s^l , one will obtain a Lind- 
blad form master equation and positivity as well as sta- 
bility of the density matrix will be granted throughout. 
In general however, this will not be the case. Here we 
will analytically consider a single qubit i/s — \ \^ ~ '^^\ 
with the simple coupling iJge — (J^ (E) B. Denoting with 
T{uj) the half-sided Fourier-transform ^TT\i of the reser- 
voir correlation function, one can write Eqn. (I12p in the 
form ps = Lps{t) and calculate the eigenvalues of the 
Liouvillian. 

If one is only interested in the subspace of diag- 
onal density matrices one finds the two correspond- 
ing eigenvalues ctq = and cti = —2X^A with 
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FIG. 5: [Color Online] Largest eigenvalue of the density ma- 
trix for one (solid) and two (dashed) qubits initialized with 
{*l Ps |j> = 1/2". calculated from Eqn. JTSJ. For largo A, not 
only positivity is violated (negative eigenvalues not shown 
- trace conservation), but the solution even becomes unsta- 
ble. The parameters in (|43|l were chosen as /3 = 1, ctJph = 1, 
cjct = 5, s = 2, such that the transition from stable to unsta- 
ble behavior occurs at Adit ~ 0.231294 for n= \. 



E. Thermalization 

In the pure dephasing limit, one observes a rapid (i.e., 
exponential) decay of the off-diagonal elements of the 
density matrix. Also, if there are no degeneracies in 
the spectrum, the rotating wave approximation (|14p pre- 
dicts a decoupled evolution of the diagonal elements of 
the density matrix in the eigenbasis of i7g, i.e., from the 
corresponding evolution equation pait) — J2j ^ij Pjj(t) 
one might expect an exponential decay into the eigen- 
vector of the matrix A which has eigenvalue (steady 
state). The corresponding subspace can still be degener- 
ate (many stationary states) or there may exist exponen- 



tially many other eigenvectors with very small eigenval- 
ues, such that thermalization does not necessarily happen 
very fast. Here we will consider some specific realizations 
of the spin-boson model ([55)1 and solve Eqn. ([H]) numer- 
ically. 

For example, if the system Hamiltonian has degenera- 
cies and the coupling to the reservoir does not lift these, 
the system may relax into states that are not even diag- 
onal in the system Hamiltonian basis. In these cases, the 
initial density matrix may determine which steady state 
is actually reached, see figure [51 As a first example we 
consider 



7 = 1, 



7i =72 



1 

2 ' 



ni = n2 = + i, 1 + i, 1 + «) : 

v2 



(70) 



and all other coefficients of Eqn. f35)) vanishing, such 
that there exists a two-fold degeneracy in the spectrum of 
iJg. In this subsection, we will consider spectral densities 
of type ([IS]) . For the low reservoir temperatures assumed 
in figure [6l thermalization corresponds to relaxation into 
the ground state and one can see that it may depend 
on the initial state of the density matrix and possibly 
lifted degeneracies whether thermalization takes place. 
In reality, the degeneracy might always be lifted by some 
imperfect Hamiltonian implementations. In addition, the 
coupling to the reservoir may be more complex than in 
([70)1 thus making the thermalized system state ([T5|) the 
only stationary state of the rotating wave master Eqn. 
(fn)l . compare also |4|]. However, it is to be expected that 
thermalization in this case will take rather long times, 
especially for large and complicated systems (see next 
subsection) . 



F. Solving Problems by Cooling 



It is known that Hamiltonians of the form ((35)) can be 
used to encode solutions to computationally hard prob- 
lems in their ground state. This is for example exploited 
in adiabatic quantum computation Q- For a system 
made of five spins, we will compare an enlarged version 
of the previous example ([70)1 



7 



5 

2 ' 



7^■ 



= 1..5 



(71) 
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FIG. 6: [Color Online] Top: Eigenvalue evolution for a 2- 
qubit system according to the secular approximation mas- 
ter Eqn. (|14p with a two- fold degeneracy for different initial 
states (solid and dashed lines) and with a near two-fold degen- 
eracy (dotted lines). Vertical dashed lines indicate the times 
at which the bottom row snapshots have been taken. 
Bottom rows: Absolute values of the matrix elements of 
the corresponding 4x4 density matrix in the ordered en- 
ergy eigenbasis (such that the top left corner corresponds to 
the ground state). Color coding ranges from blue (value 0) 
to red (value 1). In the upper row, the system is initialized 
in a non-thermal density matrix {z| pg \j) = 1/4 and relaxes 
into a the thermalized state (|15p with /3sys = /3, i.e., for the 
low reservoir temperature chosen (/3 = 10) essentially into 
the ground state. In contrast, in the middle row the initial 
state was a thermalized one with initial system temperature 



0° 

Hsy. 



1, and the system does not relax into a thermal equi- 



librium state. In the lowest row, the degeneracy was broken 
by choosing 7f = -1/2 -I- 0.01 and 7I = -1/2 - 0.01. 
Other parameters were chosen as j3 = 10, tjph = 1, uict = 5, 
s = 1, A = 0.1. 



with the ground state encoding of Exact Cover 3 - a. of the type ([35]) is given by [36|, [3 
specific NP-complete problem. 

The Exact Cover 3 problem can be introduced as 7 = to , 7^ 

follows [7]: Given a set of to constraints where each 
constraint contains the positions of three bits Cm. — 



' 2 



7»- 



(72) 



{pL^pLiPL) (with evidently 1 < pj„ < n), one is looking 
for the n-bit bit-string 5i, . . . , 6„ (with bi E {0, 1}) that 



fulfills for each constraint 6pi +bp2 



-bp3 =1 (where ' 



denotes the integer sum) . One way to encode the solution 
to this problem into the ground state of a Hamiltonian 



and all other coefhcients vanishing - the differing 
pre-factor of jfj^ in comparison to (32] results from 
the absence of double-counting in ij in Eqn. (j35p . 
In above equation, m denotes the total number of 
clauses, the number of clauses involving bit i, 
and Hij the number of clauses involving both bits 



i and j. Specifically, we will consider the 4 clauses 
Ci = (2,3,4), C2 = (1,2,5), C3 = (1,4,5), and 
C4 = (3,4,5). This implies the non- vanishing coef- 
ficients 7 = 4, 7f=7|=7| = -l, 7|=7| = -3/2, 



7ll = 7ll = 1/2, 



and 



7l5 



714 = 753 = 754 
7ll=7ll = l- 
As can be easily checked, this problem has the unique 
solution 1 10100) (with energy Eq = 0) and the first ex- 
cited states (with a six-fold degeneracy) are given by 
100101) , 110010) , 100001) , 100010) , 101001) , 101010) with 
energies Ei ^ . . . — = 1. The first excited states have 
Hamming distances (i.e., number of bit-flips necessary for 
transformation) to the solution of 2,2,3,3,4,4, respec- 
tively. This already indicates the hardness of such prob- 
lems. Simple coupling Hamiltonians such as p7p that are 
only linear in the Pauli matrices will to first order only 
yield transitions between states with Hamming distance 
1, since (a| trf \b) = if \a) and \b) have Hamming dis- 
tance larger than one. Of course, this does not completely 
prohibit transitions between states with a larger Ham- 
ming distance, but such tunneling processes will have to 
pass through energetically less favorable states and are 
therefore strongly suppressed. Accordingly, one may ex- 
pect that the process of thermalization is strongly ham- 
pered. This is also observed in figure [71 where we have 
assumed a reservoir temperature much smaller than the 
fundamental energy gap. Whereas for the simple qubit 
system (j7ip the system rapidly relaxes into the ground 
state, this is very different for the example (|7^. 



V. CONCLUSION 

We have compared different procedures of deriving 
master equations from microscopic models and have 
shown that by using a coarse-graining approach one al- 
ways obtains master equations in Lindblad form. This 
ensures for positivity and stability of the density ma- 
trix. In contrast, the usual Markovian approximation 
scheme may sometimes lead to a non-positive and even 
unstable behavior, where in the latter case there is no 
hope of approximating the exact solution. The coarse- 
grained master equations depend on a parameter - the 
coarse graining timescale r. For short coarse-graining 
times that are adaptively matched with the physical time, 
the solution Pg^r) must approximate the result of the 
Born-approximation by construction. For large coarse- 
graining times and time-independent system Hamiltoni- 
ans, we reproduce the secular approximation. For all 
intermediate coarse-graining times, a positive evolution 
of the system density matrix is ensured by the Lindblad 
form of the resulting differential equations governing the 
time evolution. For the special case of pure dephasing of 
a single qubit we reproduce the analytical result by [24| 
that p|(i) — ps{t) yields the exact solution (which is of 
course equivalent in the weak-coupling limit to the Born 
approximation) . For the simple example of exponentially 
decaying correlation functions we find by numerical simu- 
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FIG. 7: [Color Online] Top: Eigenvalue evolution for a 5- 
qubit system according to the secular approximation master 
Eqn. (|14|l for different system Hamiltonians given by Eqn. 
(|7ip (solid lines) and equation (|72p with an Exact Cover 3- 
problem (dashed lines), respectively. In both cases, the den- 
sity matrix has been initialized as {i\ pg \j) = 1/32. Whereas 
for the uncoupled spin system the system rapidly relaxes into 
the ground state, the final state for the exact cover problem 
is very different. Bottom: Absolute values of the matrix el- 
ements of the 7x7 top left sub-matrix of the full 32 x 32 
density matrix in the ordered energy eigenbasis. On the left 
(corresponding to the solid line on top), the system relaxes 
to the ground state (red spot), whereas for the exact cover 
problem (right, dashed line on top), the system state is not 
thermalized. 

Other parameters were chosen as /3 = 10, uiph = 1, i^ct ~ 5, 
s = 1, A = 0.1. 



lation a surprisingly perfect agreement between solutions 
of the integro-differential Born equation and solutions of 
the adaptive coarse-graining approach. Given an interest 
in the system density matrix at time t, we therefore pro- 
pose to match the coarse-graining time with the physical 
time T = t. In this case one has to calculate the Liou- 
villian matrix elements only once (for the desired 
time t) and then evolve the density matrix using Eqn. 
(I24p until that time t. In terms of computational com- 
plexity, this is more efficient than an evolving an integro- 
differential equation. A similar advantage is given when 
the resulting master equations are so simple such that 
an analytical solution in terms of the dampening matrix 
elements is possible. If in contrast the system is very 
complex and one is interested in the density matrix at 
all times, this advantage is destroyed. From a compu- 
tational perspective, it would therefore be interesting to 
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find the general map 



ps{t) ^ £{t)ps{t) = lim 

At^O 



psjt + At) - psjt) 
At 



which fulfills ps{t) — Ps{t)- Formally, such a map can 
be found by inserting the solution pg{t) = Vs(0) into 
the derivative in Eqn. (|73p 



Cit) 



lim — — 

At^O At 

d Ttf 

-re 

dt 



We have given various examples of simple qubit sys- 
tems coupled to a bosonic reservoir where thermalization 
which have demonstrated that thermalization of spin sys- 
tems coupled to a bosonic bath may depend on a plethora 
of factors such as the initial state, complexity of the sys- 
tem Hamiltonian and complexity of the coupling. Future 
research will consider the importance of the coarse grain- 
ing approach for different scenarios. 
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APPENDIX A: NEGLECT OF BACK-ACTION 

Throughout the paper, we make the following simpli- 
fying assumptions 

Initial factorization of the density matrix: 
By assuming 



APPENDIX B: MARKOV APPROXIMATION 

(73) 

Starting from the Born Eqn. ([8]), the Markovian ap- 
proximation in [3| is performed by replacing (in the inter- 
action picture) ps{t') — > ps{t) under the integral, substi- 
tuting T ^ t — t' and extending the integration to infin- 
ity. This is usually motivated by fastly decaying reservoir 
correlation functions By doing so, we obtain the BM 
Master Equation 

PS = -2TrB{[//sB(i),Ps(0)p°]} 

Tre { [Hs-B (t) , [Hsb (t - t) , ps (OPb] ] } 
'0{X'}, (Bl) 



(74) 



p(0)=ps(0)®p° =ps(0)®p° 



(Al) 



we implicitly demand that one can at i = prepare the 
system in a product state, which requires sufficient ex- 
perimental control. 

Born approximation: 
For a bath much larger than the system and weak cou- 
pling it is reasonable to assume that the back-action of 
the system onto the bath is small, such that the bath part 
of the density matrix is hardly changed from its initial 
value, i.e. 



Pit') = psit') • 



0{X}. 



(A2) 



Clearly, inserting this approximation in the third term in 
([7]) is consistent up to second order in A. 

Reservoir in Stationary State: 
We will assume that the reservoir is in a stationary state 
of the bath Hamiltonian, which implies [pg,i/B] = 0. 
One possibility of such a stationary state is to assume the 
thermal equilibrium state {(3~^ = ksT at temperature T) 



TrBle-z^^B} 



(A3) 



where one can now insert the decomposition ^ of the 
interaction Hamiltonian to obtain 



PS - -iXY,{BA)[AMt),Psm 



OO 

[Y, I [AB{t-T)ps{t),Aj^{t)]CAB{T)dT 

}+0{X'}, (B2) 



-h.c 



with the reservoir correlation functions ([9]). We will use 
{Ba) = throughout, since this case can always be gen- 
erated with a suitable transformation [22|. In order to 
evaluate the time-dependence of the operators in (jB2[) . it 
is useful to expand them into eigenoperators of the sys- 
tem Hamiltonian (assuming iJg to be time-independent) 



^A ' 



Aa = '^AaH = '^A\{u;) 

Aa{u;) = hE,-E^),^ \a) {a\ Aa \h) (6| , (B3) 

ah 

where the variable lo ranges over all energy differences of 
i?s, i?s I a) — Ea\a) and 5 is a Kronecker symbol. These 
operators Aa^^) have the advantageous properties 

{H^,Aa{uj)\ = -cjA^(cj), 

H^.A\{uj^ = W^(w), 

H^,A\{uj)Ab{uj'^ = 0, 
which implies in the interaction picture 



(B4) 



A^i) 



*A^(a;)=^e+-'Ai,(^). (B5) 



With inserting the half-sided Fourier transform of the 
reservoir correlation functions (jlip we obtain for the mas- 
ter equation in (jB2[) with (Ba) ~ 



Ps 



A 



{ E E e^'"''"''r^6(c.) \A,,{uj)p^{t\A\{J) 

AB ui,ui' 

+h.c.} . (B6) 
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Finally, we observe that by re-inserting the definition of 
the eigenoperators (|B3p in (|B6[) and switching back to 
the Schrodinger picture the time-dependent phase factors 
vanish and one obtains the Born-Markov master Eqn. 



APPENDIX C: SECULAR APPROXIMATION 

For large times, the terms with an oscillating prefactor 
in Eqn. (|B6p will average out and by inserting 



(CI) 



in Eqn. (|B6p and switching to the Schrodinger picture 
(compare also e.g. chapter 3.3 in Q) we obtain 



PS 



-i[Hs,ps{t)] 



Lu AB 



Ab{uj)ps,AI^{uj) 



h.c.j . 
(C2) 



The advantage of the secular approximation is that we 
can now combine the half sided Fourier transforms to 
full (even and odd) Fourier transforms (|13p of the reser- 
voir correlation functions. Inserting these definitions into 
(|C2p . one obtains a Lindblad [16] form 



PS = -■i[Hs,Ps{t)] 



2^ E E <^AB{uj)A\i^)Asiuj), ps{t) 



ui AB 

X Ab{uj)psA\{uj) - i {aI{u;)Ab{cj), ps{t)} 

(C3) 

where the positivity of the 7^ « (a;) -matrix is guaran- 
teed by Bochners theorem 0, ISg] , which states that the 
Fourier transform of a function of positive type (as are 
the reservoir correlation functions) gives rise to a positive 
definite matrix. In above equation, the second commu- 
tator corresponds to the unitary action of decoherence 
(Lamb-shift). Finally, we can insert the operator defi- 
nitions (|B3|) in jCSll to obtain Eqn. ([H]). Naturally, if 
[Hs,Ay\] — (pure dephasing), BM and BMS approxi- 
mations are equivalent. 

APPENDIX D: ADIABATIC APPROXIMATION 

With inserting the ansatz 

U{t) = ^uah(t)expi -^ J Ea{t')dt' 



ah 



X \a{t)) (6(0)1 , 



where the \a{t)) span an instantaneous basis (cho- 
sen to be complete and orthonormal) of the system 
Hilbert space defined via Hs{t) \a{t)) = Ea{t) \a{t)) in 
the evolution equation for the time evolution operator 
U{t) = —iHs{t)U{t), one obtains an equation for the ex- 
pansion cocfhcients 



Uab + Uab (a a) 



'if gca{f')dt' 



with the energy gap 

9ca{t')=E,{t')-Eait'). 

With introducing the Berry phase [s^ 
t 

la{t)^^ I {a{t')\a{t')) dt' 



(D2) 
(D3) 

(D4) 



this can also be written as 



dt 



Uabe 



-na(t)-ij gaa(t')dt' 



Ucbe 



which gives the general time evolution of the expan- 
sion coefficients Uab for any (also non-adiabatic) sys- 
tem Hamiltonian Hs{t) if one uses the initial condition 
Mab(O) — Sab- The full adiabatic approximation essen- 
tially consists in setting the right hand side of above 
equation to zero: For slowly varying system Hamilto- 
nians, the change of the eigenvectors will be negligible 
such that (a|c) « 0. Note however, that in the vicinity 
of avoided crossings, the condition of adiabaticity relates 
the maximum speed of the time evolution with the spec- 
tral properties of the system Hamiltonian 

= (D5) 

see also e.g. If (a|c) ~ 0, one obtains with the adi- 

abatic approximation u'^'^it) = SabS^'^"''-*^ which implies 
for the time evolution operator in the adiabatic limit 

t 

U{t) « Ht)){am . (D6) 

a 

Therefore, we obtain for the time-averaged operator in 
(l3l 



ah 



X- / {a{t')\Aj,\h{t')) X 



la.b(t')-^t'-J gab(r)dT 




dt' 
(D7) 



(Dl) with7,b(t') =7a(i')-7b(^') a.ndgab{T)=Ea{T)~Eb{T). 
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APPENDIX E: POSITIVITY-PRESERVING 
MASTER EQUATIONS 

Here we will show that master equations of the form 

p = -i[n{t),p{t)] 

K 

+ J2 ra^Wka(t)p(i)4W 



a/3=l 



-i{4(0/:.(t),pW} 



(El) 



generally preserve the positive semidefiniteness of an 
initial condition p(0) if the matrix Ta/iit) is positive 
semidefinite and the operator H{t) = is hermi- 

tian at all times (smoothness of all time-dependencies 
provided). With discretizing the time derivative and by 
introducing new operators 



Wi{t) 

W2{t) 



i = wl{t), 



a/3 



WK+2{t) = CK{t) 

we obtain 



K+2 



where the Wap{t) matrix is given by 



(E2) 



(E3) 



wit) 



1 


-At 


0- • -0 


-At 





0- • -0 














Atr(t) 











(E4) 



This matrix has a simple block structure and it is there- 
fore straightforward to relate the eigenvalues of w to 
those of r. In particular, one obtains the eigenvalues 
wi = 1/2 [l- x/TT4A?], W2 = 1/2 [1 + Vl + 4At2J, 
and Wi>3(t) ~ At7,;_2(^), where 7^-2 (^) are the non- 
negative eigenvalues of r{t). With diagonalizing the ma- 
trix in (jE3p via Wap(t) = 'Y^^Ua'j{t)u*p^{t)w^{t) with a 
suitable (time-dependent) unitary transformation U{t) 
we obtain from Eqn. (|E3p 



p{t + At) = Y,^7it)Wyit)Pit)W^it)^ 

7 

{t)Wc.{t). (E5) 

a 

Assuming that at time t we have a valid density matrix 
with < ps{t) < 1 we obtain by inserting the spectral 



decomposition p{t) = J2s Psi^) \'^s{t)) {<^5{t)\ 



($|p(t + Ai)|$) 



> 



> 



'^w^{t)ps{t) 

7(5 



- VI + 4At2 

J2ps{t)\{^\m\^s{t)) 



At-+0 

> 0. 



[<^>\Wi\<S>s{t)) 

(E6) 



such that in the limit At ^ (which defines the origi- 
nal differential Eqn. (jEip ). the smallest eigenvalue of the 
density matrix at time t + At approaches zero faster than 
the discretization width At. Therefore, in this limit the 
matrix w(t) becomes positive semidefinite and the differ- 
ential Eqn. (jE3[) becomes a positivity-preserving map. 



APPENDIX F: SING DISTRIBUTION 

For discrete a, b and continuous oj we would like to 

analyze 

/(w,a, 6) = lim r sine {uj + a)— sine (u + b) — 
r^oo L 2J L 2. 

= ,^4sin[(. + a)^]sin[(^ + 5)^] 

r^oo {to + a){uj + b)T ^ ' 

One can consider the case a ^ 6 by a partial fraction 
expansion where without loss of generality we find for 
the first term (due to the symmetry the second term can 
be treated in a completely analogous way) 



-t-00 



COS 



sin[(tj + a)§] sin[(w + 6)§] 

g{uj) ^ — \ ^dw 

T\LO + a) 



(h 7 ( ^ sin^ [(a;-fa)f] 

2 J J T{uj + a) 



-— sm 
2 



+ 00 

Ti [ sin [(cj + a)T] 
(b — a)— I q(ijj) r — dio . 

— CO 



where we have inserted {uj + b) ~ {uj + a) + {b— a) to use 
trigonometric relations for sin(a + /?). With a suitable 
transformation, this becomes 



cos [(6 — a)§] 



3in [(6 — a)-, 
2^ 



+00 



X \ sm (x 2) , 
g I a] dx 

T I X 



sin(a;) 



dx j(F2) 
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which yields Eqn. (|27|) . 



hm r sinc^ [Qr] = n5{n) (F3) 

r— *oo 

one can consider the general case via 
f{uj^a^h) X 5ab lim t sinc^ 

r—>oo 

= (^) , (F4) 



where for large r, the dominant integral contribution 
evaluates the (smooth) function near g{—a) such 

that this value can be taken out of the integral and we 
obtain lim la^b ~ 0. Using the representation 

r — >oo 
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